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Abstract 

Heterogeneous anisotropic diffusion problems arise in the various areas of science and en- 
gineering including plasma physics, petroleum engineering, and image processing. Standard 
numerical methods can produce spurious oscillations when they are used to solve those prob- 
lems. A common approach to avoid this difficulty is to design a proper numerical scheme and/or 
a proper mesh so that the numerical solution validates the discrete counterpart (DMP) of the 
maximum principle satisfied by the continuous solution. A well known mesh condition for the 
DMP satisfaction by the linear finite element solution of isotropic diffusion problems is the non- 
obtuse angle condition that requires the dihedral angles of mesh elements to be non-obtuse. In 
this paper, a generalization of the condition, the so-called anisotropic non-obtuse angle condi- 
tion, is developed for the finite element solution of heterogeneous anisotropic diffusion problems. 
The new condition is essentially the same as the existing one except that the dihedral angles are 
now measured in a metric depending on the diffusion matrix of the underlying problem. Several 
variants of the new condition are obtained. Based on one of them, two metric tensors for use in 
anisotropic mesh generation are developed to account for DMP satisfaction and the combination 
of DMP satisfaction and mesh adaptivity. Numerical examples are given to demonstrate the 
features of the linear finite element method for anisotropic meshes generated with the metric 
tensors. 
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1 Introduction 

We are concerned with the numerical solution of the diffusion equation 

-V-(DVu) = /, in Q (1) 
subject to the Dirichlet boundary condition 

u = g, on dft (2) 
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where flcR (d=l,2, or 3) is the physical domain, / and g are given functions, and D = D(ac) 
is the diffusion matrix assumed to be symmetric and strictly positive definite on fi. The boundary 
value problem (BVP) ([!]) and ^ becomes a heterogeneous anisotropic diffusion problem when D 
changes from place to place (heterogeneous) and its eigenvalues are not all equal (anisotropic) at 
least on a portion of f2. When a standard numerical method, such as a finite element, a finite 
difference, or a finite volume method, is used to solve this problem, spurious oscillations can occur 
in the computed solution. A challenge is then to design a proper numerical scheme and/or a proper 
mesh so that the computed solution is free of spurious oscillations. In some applications such 
as plasma physics, it is further desired that the mesh be aligned with the fast diffusion direction 
so that no excessive numerical dissipation is introduced in slow diffusion directions. Moreover, 
mesh adaptation is often necessary for improving computational efficiency and accuracy when the 
physical solution and/or the diffusion matrix have sharp jumps. 

Anisotropic diffusion problems arise in the various areas of science and engineering including 
plasma physics in fusion experiments and astrophysics [251 IZE E?l EH EU [63] , petroleum reservoir 
simulation [H [21 [TH [22J [53], and image processing O [131 E3 EH EH EH]- In plasma physics, 
magnetized plasmas are constrained to move primarily along magnetic field lines. Their heat 
conductivity in the direction parallel to the magnetic field is much higher than those perpendicular 
to it, and the ratio of the conduction coefficients can easily exceed 10 10 in fusion experiments. 
The numerical simulation of the heat conduction of plasmas must not only produce a physically 
meaningful temperature distribution but also avoid excessive numerical dissipation in the directions 
perpendicular to the magnetic field. In petroleum engineering, fluids such as water, crude oil, and 
natural gas are stored in reservoir rocks filled by interconnected networks of pores. The diffusion 
and flow of those fluids depend crucially on the rocks' permeability which changes with location 
and flow direction and has much large values in horizontal directions than in the vertical direction. 
Finally, PDE-based anisotropic diffusion filters have been successfully used for shape recognition 
and edge detection in image processing. 

The BVP and ^ is a representative example of anisotropic diffusion problems arising in 
those areas. As typical for diffusion problems, it satisfies the maximum principle 

max u(x) < maxjO, max q(s)} (3) 
xtnudn v ' ~ L sean 

provided that f(x) < holds for all x G £1. The BVP has been studied extensively in the past, 
and a major effort has been made to avoid spurious oscillations in the numerical solution. A 
common strategy is to develop numerical schemes satisfying the discrete counterpart of ^ - the 
so-called discrete maximum principle (DMP), which are known to produce numerical solutions free 
of spurious oscillations [TH EH]- The studies can be traced back to early works by Varga [68J, 
Ciarlet |14] , Ciarlet and Raviart [15] , and Stoyan [641 [65] where a number of sufficient conditions in 
a general and abstract setting are obtained for a class of linear elliptic partial differential equations 
(PDEs). For example, denote by Au = f the linear algebraic system resulting from the application 
of a numerical scheme to a linear elliptic PDE supplemented with a Dirichlet boundary condition, 
where A is the n x n stiffness matrix, u is the unknown vector, and / the right-hand-side vector. 
Then, a sufficient condition is given as follows. 
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Lemma 1.1 (I65\j ) If the stiffness matrix A satisfies 

(a) that A is monotone with A^ being either nonsingular, or singular and irreducible; and (4) 

(b) that A^e^ >0, (5) 

then the numerical scheme satisfies DMP. 

Here, matrix A is said to be monotone if A is nonsingular and A~ l > (i.e., all entries of A' 1 
are non-negative), and ' and are defined as 



a { - ] 



an, fori=j 

aij, for i / j, < , e (n) 
0, for i / j, aij > 



(6) 



Note that condition Q is equivalent to that A^ has nonnegative row sums. Moreover, A = A^ 
and the condition Q holds when A is an M- matrix |67| . From Lemma we have the following 
lemma. 

Lemma 1.2 If the stiffness matrix A is an M-matrix and has nonnegative row sums, then the 
numerical scheme satisfies DMP. 

Numerical schemes satisfying DMP have been developed along the line of those sufficient con- 
ditions by either designing a proper discretization for the underlying PDE or employing a suitable 
mesh. To date most success has been made for the isotropic diffusion case where D is in the 
scalar matrix form, D = a(x)I, with a(x) being a scalar function; e.g., see [H [TOj [151 EH EH [381 
|4"4"1 l4"9l [66l 171] . In particular, it is shown in [9l [15] that the linear finite element method (FEM) 
satisfies DMP when the mesh is simplicial and satisfies the so-called non-obtuse angle condition 
requiring that the dihedral angles of all mesh elements be non-obtuse. In two dimensions this 
condition can be replaced by a weaker condition (the Delaunay condition) that the sum of any 
pair of angles opposite a common edge is less than or equal to ir [49, 66J. Similar mesh conditions 
are developed in |36l 137] IBB] 144] for elliptic problems with a nonlinear diffusion coefficient in the 
form D = a{x, u, Vu)7 and with mixed boundary conditions. Burman and Ern [10] propose a non- 
linear stabilized Galerkin approximation for the Laplace operator and prove that it satisfies DMP 
on arbitrary meshes and for arbitrary space dimension without resorting to the non-obtuse angle 
condition. 

On the other hand, the anisotropic diffusion case is more difficult and only limited success 
has been made [HI [HI [251 E3 [23 H31 HSl HH SSI EQl EH E21 [531 EI]. For example, Draganescu 
et al. [TH] show that the non-obtuse angle condition fails to guarantee DMP satisfaction in the 
anisotropic diffusion case. The techniques proposed by Liska and Shashkov [S2] and Kuzmin et al. 
|43] to locally modify (or repair) the underlying numerical scheme, by Sharma and Hammett [61] to 
employ slope limiters in the discretization of the PDE, by Mlacnik and Durlofsky [53] to optimize 
the mesh for a multipoint flux approximation (MPFA) finite volume method (e.g., see [HE] for the 
method), and by Li et al. [50] to optimize a triangular mesh for the finite element solution, help 
reduce spurious oscillations. A nonlinear, first order finite volume method developed by Le Potier 
|4T] and further improved by Lipnikov et al. [51] gives rise to a stiffness M-matrix on arbitrary 
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meshes when applied to parabolic PDEs but fails to satisfy DMP when applied to steady-state 
elliptic problems. A first order finite difference method having similar features is proposed by Le 
Potier 03]. 

In this paper we study the linear finite element solution of BVP ([I]) and ^ with a general 
diffusion matrix D = D(x). The objective is threefold. The first is to develop a generalization of 
the well known non-obtuse angle condition, the so-called anisotropic non-obtuse angle condition 
(cf. equation (24)), so that the linear FEM satisfies DMP when the mesh is simplicial and satisfies 
this condition. The condition requires that the dihedral angles of all mesh elements, measured in a 
metric depending on D, be non-obtuse. It reduces to the non-obtuse angle condition for isotropic 
diffusion matrices. It also reproduces several existing mesh conditions for homogeneous anisotropic 
media for which D is a full, constant matrix (see Remark 2.2). The second objective is to derive 
a metric tensor for use in mesh generation based on the anisotropic non-obtuse angle condition. 
This is done by adopting the so-called M-uniform mesh approach [31] where an anisotropic mesh is 
generated as an M-uniform mesh or a uniform mesh in the metric specified by a tensor. M-uniform 
meshes generated with the metric tensor satisfy the anisotropic non-obtuse angle condition and are 
aligned with the diffusion matrix O (cf. jQ. The final objective is to combine both mesh adaptivity 
and DMP satisfaction in the numerical solution of anisotropic diffusion problems. An optimal met- 
ric tensor (see (55)) accounting for both considerations is obtained by minimizing an interpolation 
error bound, and advantages of using adaptive, DMP-bound meshes are demonstrated in numer- 
ical examples. To the authors' best knowledge, this is the first effort that mesh adaptivity and 
DMP satisfaction are considered simultaneously in the numerical solution of anisotropic diffusion 
problems. 

The outline of this paper is as follows. In section [2j the linear finite element solution of and 
((2]) is described and the anisotropic non-obtuse angle condition and several variants are derived. 
Section[3]is devoted to the derivation of the metric tensor based on the anisotropic non-obtuse angle 
condition. In section [4j the combination of mesh adaptation and DMP satisfaction is addressed, 
and an optimal metric tensor is obtained by minimizing an interpolation error bound. Numerical 
examples are presented in section [5} Finally, section [6] contains conclusions and comments. 



2 Anisotropic non-obtuse angle conditions for linear finite element 
approximation 

Consider the linear finite element solution of BVP ([I]) and (|2j). Assume that Q is a connected 
polygon or polyhedron and an affine family of simplicial triangulations {Th} is given thereon. Let 

U g = {veH 1 (n)\v\ dn = g}. 

Denote by Ug C U g the linear finite element space associated with mesh Th- Then a linear finite 
element solution u h G U g to BVP |l| and Q is defined by 

[ (Vv h ) T BVu h dx = ( fv h dx, V/ G [/ h (7) 
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where Uq = Ug with g = 0. This equation can be rewritten as 

E [ (V/) r DV^= E [ f vh dx, \fv h eU^. (8) 

Generally speaking, the integrals in ^ cannot be carried out analytically, and numerical quadrature 
is needed. We assume that a quadrature rule has been chosen on the reference element K for this 
purpose, 

„ mm 

" " (9) 



n III, lit 

k=l k=l 



where w^s are the weights and b k s the quadrature nodes. A 2D example of such quadrature rules 
is given by Wk = 3 (k = 1,2,3) and the barycentric coordinates of the nodes (g, g, |), (4, |, g), and 
(|, g, |); and a 3D example is W{ = | (i = 1,2,3,4) and the barycentric coordinates of the nodes 
(a, a, a, 1 — 3a), (a, a, 1 — 3a, a), (a, 1 — 3a, a, a), and (1 — 3a, a, a, a) with a = 5 ~^^ ; e.g., see [2~T] . 

Let Fx be the afhne mapping from K to K such that -ftT = Fk{K), and denote 6^ = Fif (£>&), 
k = 1, ■ ■ ■ ,m. Upon applying Q to the integrals in Q and changing variables, the finite element 
approximation problem becomes seeking u h E Ug such that 

m m 

E |K|^<MV W V) T ©(6f)VuV= E l#lE^/(^)^)> V ^ €[/ o (10) 
iCeT^ fc=l KeT h k=i 

where Vv h \x and S7u h \K denote the restriction of \7v h and \7u h on iT, respectively. Note that we 
have used in (10) the fact that Vv h \x and Vu h \x are constant. Letting 

m 

ox = E^ D (^)' ( n ) 
k=i 



we can rewrite (10) into 



We now express (12) in a matrix form. Denote the numbers of the elements, vertices, and 
interior vertices of Th by N, N v , and N V i, respectively. Assume that the vertices are ordered in 
such a way that the first N V i vertices are the interior vertices. Then Uq and u h can be expressed as 

Uq = spamj>i, • • • , <f) Nv .} (13) 

and 

u h =Y /Uj <j )j + e u ^ ( 14 ) 

j=l j=N vi +l 

where <pj is the linear basis function associated with the j'-th vertex, aj. Note that the boundary 
condition ([2]) can be approximated by 

Uj = Qj = g(aj), j = N vi + 1, N v . (15) 
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Substituting (14) into and taking v h = (pi (i = l,...,N v i) 
equations with (15), we obtain the linear algebraic system 

Au = f, 

where 



in (12) and combining the resulting 

(16) 



A 



An A X2 
I 



(17) 



/ is the identity matrix of size (N v — N v i), and 

u = (ui, ...,u Nm ,u Nm+ i, ...,u Nv ) T , 

f = (/l) ■••) fN vi ,9N vi +l: ■■■,9N V ) T - 

The entries of the stiffness matrix A and the right-hand-side vector / are given by 

aij = \K\ (V4>i\ K ) T O k V(j)j\ K , i = 1, .-,N vi , j = 1, ...,N V , 
KeT h 



(18) 
(19) 



fi= E I^E^)^^)' i = l,...,N vi . 

KeT h k=i 

We recall that (16) and (17) have been obtained under the Dirichlet boundary condition Q. It is 
not difficult to show that a linear system in the same form can be obtained for mixed boundary 
conditions provided that To ^ 0, with being the part of the boundary where the Dirichlet 
condition is imposed. Therefore, the mesh conditions developed below also work for mixed boundary 
conditions with ^ 0. 

We now study under what mesh conditions the scheme (16) satisfies DMP. Our basic tool is 
Lemma 1.2, i.e., we show that A is an M-matrix and has non-negative row sums when the mesh 
satisfies the condition (24) below. To this end, we first introduce some notation. Denote the vertices 
of K by af - , a^, • • • , af +1 . The edge matrix of K is defined as 



Ek 



a 



K 



a l , a 3 



a 



K 

1 ! 



) a d+l 



a 



A'l 



From the definition of simplices, Ek is nonsingular |62j . Then, a set of q- vectors (cf. Fig. [TJ can 
be defined as 



E 



K ' 



d+l 
i=2 



(20) 



This set of vectors has the following properties, 
(i) By definition, it follows that 



K 



if ■ (af 



af ) = Sij, 



a 



,d + l; j = !,-■■, d+l 



(21) 



where Sij is the Kronecker delta function. 



(ii) Denote by Sf the face opposite to vertex af (i.e., the face not having as a vertex). Then 
(21) implies that qf is the inward normal to the face Sf; see Fig. [lj 
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(iii) The dihedral angle, a^, between any two faces Sf and Sf (i ^ j) is defined as the supplement 
of the angle between the inward normals to the faces. It can be calculated by 

qf . q K 

cos(aij) = - ' 3 K i / j. (22) 
\\Qi II Wlf II 

(iv) It is known [U E] that, for any vertex of K with the global and local indices i and ix, 
respectively, there holds 

V<t>iW = q? K . (23) 




Figure 1: A sketch of the q vectors for an arbitrary element. The angles sharing the edge connecting 
vertices a\ and a 2 are a and /3. 

The main result of this section is stated in the following theorem. 



Theorem 2.1 If the mesh satisfies the anisotropic non-obtuse angle condition 
(qf) T B K qf <0, Vi^j, i,j = l,2,...,d+l, VK G T ft 
i/ien i/ie linear finite element scheme (12) for solving BVP |I]) and |Ip satisfies DMP. 



(24) 



Proof. We prove this theorem using Lemma 1.2 That is, we show that the stiffness matrix A 



has non- negative row sums and is an M- matrix when the mesh satisfies condition ( 24 ) . 



(i) We first show that A has non-negative row sums. From (17) we only need to show X^?=i a ij — 
for i = 1, ...,N v i. From (18) we have 



N v N v _ 

Yl ai i = Yl \ K \ (y<Pi\ K ) T B K 

j=l j=i K &T h 

I N v 

K£T h \j=l 

= o, 

where we have used the fact that Y^,f=i (fijfa) = 1 f° r an y x E K. 



K 



(25) 
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(ii) Next we show that 



dij < 0, 
an > 0, 



V i + j, i,j = 1, ...,N V 
Vi = l,...,N v . 



(26) 
(27) 



Let uji (or ojj) be the patch of the elements containing otj (or a,j) as a vertex. Notice that S7(pi\K = 
when K ^ uji. Denote the local indices of vertices ctj and a, on K by %k and j^-, respectively. Then 
from (18), (23), and (24), we have, for i ^ j, i 



l,...,N vi , j 
\K\ (V(t>i\ K ) T 



< 



0. 



1,...,N V , 

V4>j\ K 



hi q 



K 
3K 



N vi + 1,---N V , j 



(28) 
(29) 

1,...,A^. Hence, the 



From (17) it is obvious that = for i ^ j, : 
off-diagonal entries of A are non-positive. 

The inequality (27) follows immediately from (17), (18), and the positive defmiteness of I$k- 
(iii) We now show that An defined in (17) is an Af-matrix. Notice that the non-negativeness 
of the row sums of A and the properties (26) and (27) imply that An is diagonally dominant. In 
theory, we can show that An is an M-matrix by proving it is irreducible [67]. However, we will 
need to assume that any pair of interior vertices is connected at least by an interior edge path [18] . 
To avoid this additional restriction on the mesh, we instead opt to show An is symmetric and 
positive definite, which together with (26) and (27) implies that An is an M-matrix [67] . 

From (18) it is obvious that An is symmetric. It suffices to show An is positive definite. From 
the strictly positive definiteness of the diffusion matrix D, there exists a positive constant /3 such 
that 

B K > (31, VK G T h . 



For any vector v 



N - 

(vi, VN vi ) , we define v = ^ Vi<j>i G Uq . From the definition of An and the 



fact that \7v \k is constant on K, we have 



v T Anv 



\K\ (Vv h \ K f B K Vv h \ K 

KeT h 



k) T Vv h \ K 



(Vv h ) T Vv h dx 



KeT h 



K 



= p j (Vv h ) T Vv h dx 

> pc p 

where in the last step we have used Poincare's inequality and C p > is the associated constant. For 



\v h \ 2 dx, 
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any nonzero vector v, v 



N vi 
i=i 



^ and is piecewise linear and continuous on fi. Consequently, 



v T A u v > (3C P [ \v h \ 2 dx > 0, V« ^ 
which implies that An is positive definite. Hence, An is an M-matrix. 



(iv) From (17) it is easy to verify that the inverse of A is given by 



A- 



A-^ 1 -A^An 
I 



Then (26) and the fact A^ > imply that A 1 > and therefore A is an M-matrix. 



We have shown above that A is an M-matrix and has non- negative row sums. By Lemma 1.2 

□ 



we conclude that the linear FEM satisfies DMP when the simplicial mesh satisfies (24). 



Remark 2.1 For the isotropic case where D = a(x)I for some scalar function a(x), condition 



( 24 ) reduces to the well known non-obtuse angle condition IS] 



qf-qf<0, Vi^j,VK£T h , 



(30) 



which requires the dihedral angles otij (cf. (22)) of all mesh elements be non-obtuse. Thus, condition 



(24) is a generalization of the non-obtuse angle condition. An alternative interpretation of (24) is 



that the dihedral angles of element K, measured in the Riemannian metric (piecewise constant), 
are non-obtuse. □ 



Remark 2.2 It is interesting to point out that an explicit mesh condition similar to (24) is 



obtained by Eigestad et al. [20] for a multipoint flux approximation (MPFA) finite volume method 
on triangular meshes for anisotropic homogeneous media (i.e., O is constant). Moreover, (24) 
reduces to a mesh condition obtained by Li et al. |50| for a similar situation with constant O and 
triangular meshes. To see this, let the eigen-decomposition of the constant diffusion matrix D be 



cos ( 
sin( 



- 8111 1 

cos l 



h 
k 2 



cos ( 
- sin( 



sm( 
cos ( 



(31) 



For an arbitrary triangular element K, denote the angles sharing the edge connecting vertices a\ 
and ci2 by a and /3; see Fig. [T] Then, a mesh condition of [50] is given by 



-k\ sin f3 sin a + k 2 cos (3 cos a < 0, 
-&2 cos j3 < 0, 
-&2 cos a < 0, 



(32) 



provided that the edge connecting a\ and a 2 is parallel to the primary diffusion direction (cos 9, sin 9) T 
(the eigenvector corresponding to the first eigenvalue of D, k\). We now show that (24) reduces 



to (32) for the current situation. Without loss of generality we assume that the primary diffusion 
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direction and the edge connecting a\ and a 2 are in the direction of the x-axis; cf. Fig. [IJ (In this 
case we have 9 = 0.) It is not difficult to obtain 



Qi = ci 



- sin/3 
cos (3 



<?2 = c 2 



sin a 

cos a 



<?3 = c 3 




1 



where ci, c 2 , and C3 are positive constants. From these and (31), (24) reduces to 

qJH>Kq 2 = Qi^Qz = c \ c 2{—k\ sin a sin /3 + &2 cos a cos /3) < 0, 
qjB K q 3 = qfOq 3 = ac 3 (-k 2 cos f3) < 0, 
q^KQz = Q2 n <l3 = c 2 c 3 (-k 2 cosa) < 0, 



which gives ( 32 ) . 



It is often more convenient to express the anisotropic non-obtuse angle condition ( 24 ) in terms 



of mapping Fk from K to K. Denote the Jacobian matrix of Fk by F' K . We define the vectors q k , 
k = 1, d + 1 for the reference element K as in (20). The chain rule of differentiation implies 



= (F' K )- T V^ h 



where = 4>i(Fk(^)). From (23), we have 



(F' K )- T qi . 



Inserting this into (24) we obtain the following theorem. 



Theorem 2.2 If the mesh satisfies 



qf (F' K )- 1 B K (F' K )- T q ] < 0, Vi^j, i, j = 1, a! + 1, VXeT A 



(33) 



i/ien i/ie linear finite element scheme (12) for solving BVP \1V and |5p satisfies DMP. 



Corollary 2.1 Suppose that the reference element K is taken as a simplex with non-obtuse 
dihedral angles. If the mesh satisfies 

- T = C K I, VK G T h 



{F' k )- 1 ®k{F' k )- 



ih (34) 
where Ck is a positive constant on K and I is the d x d identity matrix, then the linear finite 



element scheme (12) for solving BVP |7l) and |M) satisfies DMP. 



Proof. Since K is a simplex with non-obtuse dihedral angles, we have 

i,j = l,...,d+ 1. 



qf qj < 0, * ^ j, 



From this it is easy to see that (34) is sufficient for (33) to hold. 



In the next two sections mesh condition (34) will be used to develop metric tensors accounting 



for DMP satisfaction and mesh adaptivity. These metric tensors are needed in anisotropic mesh 



generation. It is emphasized that (34), as well as mesh conditions (24) and (33), can also be used 



more directly via direct minimization |5Q|, 153] or variational formulation |30j for optimizing the 
current mesh to improve DMP satisfaction. 
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3 Anisotropic mesh generation: metric tensor based on DMP sat- 
isfaction 



In this section we develop a metric tensor for use in anisotropic mesh generation based on mesh 
condition (34). To this end, we adopt the so-called M-uniform mesh approach [311 [32] where an 
anisotropic mesh is viewed as an M-uniform mesh or a uniform one in the metric specified by a tensor 
M = M(x). The tensor, chosen to be symmetric and positive definite, provides the information on 
the size, shape, and orientation of mesh elements over Q necessary for the actual implementation 
of mesh generation. Various formulations of the metric tensor have been developed in the past for 
anisotropic mesh adaptation; e.g., see [H HH [231 ED S3] ■ Once a metric tensor has been determined, 
the corresponding anisotropic meshes can be generated using a variety of techniques including 
blue refinement [121 05], directional refinement [581 [59], Delaunay-type triangulation [3] IS"! [TT] 156]. 
front advancing [21], bubble packing [72] , local refinement and modification [31 [H [TTJ, [23 [60] , and 
variational mesh generation [71 [191 E01 EHl HD US EQ] ■ In this paper we restrict our attention to the 
determination of a metric tensor for DMP satisfaction, and refer the interested reader to the above 
mentioned references for meshing strategies. 

It is shown in |32| that when the reference element K is taken to be equilateral and unitary in 
volume, a simplicial M-uniform mesh Th for a given M = M(x) satisfies 

crh 



PK\K\ 



N ' 



\/K G Th 



tr((F' K ) T M K F' K ) = & e t((F' K ) T M K F' K y , VKeT h 



(35) 
(36) 



where N is the number of mesh elements, Fk is the affine mapping from K to K, F' K is the Jacobian 
matrix of Fk, and 

1 



M K = TFT f M(x)dx, PK = v / det(M x ) 
\ K \ Jk 



KeT h 



(37) 



Condition (35), referred to as the equidistribution condition, determines the size of K from pk- 
The larger pK is, the smaller \K\ is. On the other hand, (36), called the alignment condition, 
characterizes the shape and orientation of K in the sense that the principal axes of the circumscribed 
ellipsoid of K are parallel to the eigenvectors of Mr while their lengths are reciprocally proportional 
to the square roots of the respective eigenvalues [32J . 

To determine M from mesh condition (34), we first notice that the left and right sides of 
(36) represents the arithmetic and geometric means of the eigenvalues of matrix (F' K ) T Mj<F' K , 
respectively. From the arithmetic- mean geometric- mean inequality, (36) implies that all of the 
eigenvalues are equal to each other. In other words, {F' K ) T MkF' k is a scalar matrix, i.e., 



{F' K ) T M K F' K = C K I or (Fx)* 1 Mft 1 (F' K j 



-T 



Ck 1 ! 



(38) 



for some constant Ck- A direct comparison of (38) with (34) suggests that the metric tensor M 
be chosen in the form 

M DM p,K = e K ^ K \ VK£T h (39) 

where 9 = 9k > is an arbitrary piecewise constant function. Thus, any M-uniform mesh associ- 
ated with a metric tensor in the form (39) satisfies condition (34). The following theorem follows 



from Corollary 2.1 
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Theorem 3.1 Suppose that the reference element K is taken to be equilateral and unitary in 



volume. For an M -uniform mesh associated with any metric tensor in the form (39), the linear 
finite element scheme (12) for solving BVP Hy and M) satisfies DMP. 



Remark 3.1 Since an M-uniform mesh is aligned with the metric tensor M as characterized 



by the alignment condition (36), we can conclude that when M is chosen in the form (39), a 



corresponding M-uniform mesh is aligned with the diffusion matrix D in the sense that the principal 
axes of the circumscribed ellipsoid of element K are parallel to the eigenvectors of Hk while their 
lengths are proportional to the square roots of the respective eigenvalues. As a consequence, the 
length of K is greater in a faster diffusion direction and smaller in a slower diffusion direction. A 
small length scale of mesh elements in slow diffusion directions helps reduce numerical dissipation 
in those directions. □ 



Remark 3.2 Note that 9 = 9 k in (39) is arbitrary. Thus, in addition to satisfying DMP, there 



is a degree of freedom for the mesh to account for other considerations. In the next section we shall 
consider mesh adaptation and choose 9k to minimize a certain error bound. □ 



4 Metric tensors based on DMP satisfaction and mesh adaptivity 

In this section we develop a metric tensor taking both the satisfaction of DMP and mesh adaptivity 



into consideration. The metric tensor takes the form (39), with the scalar function 9 = 9k being 



determined to minimize an interpolation error bound. For simplicity, we consider here an error 
bound for linear Lagrange interpolation. Other interpolation error bounds (e.g., see |32j) can be 
considered without major modification. 

Lemma 4.1 (132f ) Let K C M rf be a simplicial element and Tlh be the linear Lagrange interpo- 
lation operator. Then, 



U h v\ m{K) < C\\(F' K )- l \\ \J [tr((F' K ) T \H(v)\F' K )} 2 dx 



G H 2 (K) 



(40) 



where || • || denotes the I2 matrix norm, H(v) is the Hessian of v, and \H(v)\ = \J H(v) 2 . 
Lemma 4.2 For any given d x d symmetric matrix S, there holds that 

\tr(A T SA)\ < tr(A T A) ||S||, G R dxd . (41) 
// S is further positive definite, then 

\\S\l~ 1 tr(A T SA) < tr(A T A) < tr(A T SA) \\S~\ (42) 
Proof. Denote the eigen-decomposition of S by 

S = QT,Q T , 

where Q is an orthogonal matrix, S = diag(Ai, A^), and Aj, i = 1, d are the eigenvalues of S. 
Write 

A T Q = [ Vl ,...,v d ). 
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Then 



It follows that 



A T SA = (A T Q)Z(Q T A) = [ Vl , v d ]E[v 1} ...,v d f = ^ X^v] 

i 

\t?{A T SA)\ = |^A;tr(^f)| 

i 

= |^Ai||^i|| 2 | 

i 

■max 

i 

= tr(A T A) \\S\\, 



which gives (41). Inequality (42) follows from (41) and that 

tr{A T A) = tr(A T S^S- 1 S^A) < tr(A T SA) \\S' l \ 



(43) 



The scalar function 6 = Ok in ( 39 ) is determined based on interpolation error bound ( 40 ) . From 
the definition of the Frobenius matrix norm, we have 



A\\ < \\A\\ F = Jtr(^ T ^) = Jtr(^ T ), VA G 



odxd 



Using this, taking squares of both sides of (40), and summing the result over all elements of 7~h, we 
have 



\u-U h u\ 2 H1{n) = \ u -n h u\ 2 Hl{K) 
KeT h 

<CY1 WiFkr'f [ [tr((F' K ) T \H(u)\F' K )} 2 dx 
Ken Jk 

<CY. \\( F Kr l \\l [ [tT((F' K f\H(u)\F K )] 2 dx 



Ken 



Ken 

From Lemma I4J2] it follows that 



H(F'K)-\F' K r T )} [ [tr ((F' K f\H(u)\F' K )] 2 dx. 

JK 



\u - n h «|^i (n) 

[tr((fJ c )- 1 D J r(i^)- T )]-[|D*l|- / [tr((F^ T B K \F' K ))] 2 \\B K \H(u)\fdx 
Ken Jk 



CY,\K\- [tT((F' K )- l O K (F' K )- T )] • [tr((^W(^))] : 



Ken 



^II-T^T I \pK\H(u)\fdx. 
\ K \ Jk 



(44) 
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Consider an M- uniform mesh Th corresponding to a metric tensor Mk in the form (39). Then, 



alignment condition (36) reduces to 

1 



d 



tr {(F' K fD- K l F' K ) = det ((i^) T D^i^) a 



(45) 



From the arithmetic- mean geometric- mean inequality, (45) implies that all of the eigenvalues of 
matrix (F , k ) t H>k F' K are equal to each other. As a consequence, all of the eigenvalues of the 



inverse of (F' K ) T H K F' K are equal to each other, which in turn implies 



1 



tr {{F' K )- l H K {F' K r T ) = det {{F' K )" 1 H k{F' k )~ T ) 3 . 



(46) 



Inserting (45) and (46) into (44) and noticing 



det ((F^fO^F^) = |K| 2 det (D^)" 1 , det {{F' K )- X H K {F' K 



/ \-T\ 



we have 



\u-U h u\ 2 Hl{n) < C |-^|~det(I 



K) 



Rewrite this bound as 



where 



K 



\K\ 



K|- 2 det (B K ) 
H(u)\\\ 2 dx. 



K 



u ■ 



n hU \ 2 Hl(u) <c \k\—b k , 

KeT h 



B K = det (B K y* 



ft — 1 1 



1 



K 



D K \H(u)\\\ 2 dx. 



(47) 



(48) 



(49) 



\K\ 

Notice that f K \\U)k\H(u)\\\ dx and therefore Bk can vanish locally. To ensure the positive defi- 
niteness of the metric tensor to be defined, we regularize the above bound with a parameter ah > 
as 

,o ^ — •>, . , d+2 ^ — -v , d+1 

\u-U h u\ 2 Hl(n) <C 2^ \K\ d [a h + B K ] = Ca h \ K \ d 
KeT h Ker h 
From Holder's inequality, we have 



1 + — B K 

ah 



(50) 



El* 

KeT h 



d+2 
d 



1 + —B K 

ah 




l + —B K 

ah 



d ' 
d+2 



> N~d ^ \K\ 



1 + —B K 

ah 



d+2 
d 



d+2 



d+2 
d 



(51) 



with equality in the last step if and only if 



\K\ 



l + —B K 

ah 



d 

d+2 



constant, MK 6 Th- 



(52) 



A direct comparison of this with equidistribution condition ( 35 ) suggests that the optimal pk be 
defined as 



PK 



1 + — B K 

ah 



d 
d+2 



(53) 
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From the relation px = \J det(Mx), we find the optimal Ok and Mk as 



9 K = p K det (Da-) a 



1 + — 5 X 

"ft 



2 

d+2 i 

det (D^) d 



M 



DMP+ adap,K 



l + —B K 



2 

d+2 



det 



»K) d ^ K 



(54) 



(55) 



where Bk is defined in (49). With the so-defined metric tensor, the error bound can be obtained 



from (50) and (51) for a corresponding M-uniform mesh as 



d+2 
2d 



(56) 



To complete the definition, we need to determine the regularization parameter ah- We follow 
31] to define ah such that 

<Th= PK\K\<2\n\, (57) 



KeT h 



with which roughly 50% of the mesh points are concentrated in regions of large px- From (53) and 
Jensen's inequality, we have 



07, 



= Em 

KeT h 

< £1*1 

KeT h 



d 
d+2 



1 + — B K 

ah 



l + a h d + 2 B K + 2 



>d+2 

>K ■ 



= N + «/ +2 E \ K \ B i 
Kar h 

By requiring the above bound to be less than or equal to 2|f2|, we obtain 



(58) 



an 



d 

I d+2 

>K 



d+2 
d 



(59) 



KeT h 



Combining (56) with (57) and (59) and summarizing the above derivation, we have the following 
theorem. 

Theorem 4.1 Suppose that the reference element K is chosen to be equilateral and unitary in 



volume. For any M-uniform simplicial mesh corresponding to the metric tensor (55), the linear 



finite element scheme (12) for solving BVP and IM) satisfies DMP and the interpolation error 



for the exact solution u is bounded by 



d+2 
2d 



|«-n^| H i (n) < cn d 



L \ k \ b k 2 

J<eT h 



(60) 



where Bk is defined in (49). 
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It is remarked that the metric tensor (55) (cf. (49)) depends on the second derivatives of the 



exact solution u which is what we are seeking/approximating. In actual computation, the second 
derivatives are replaced with approximations obtained with a Hessian recovery technique such as 
the one of using piecewise quadratic polynomials fitting in least-squares sense to nodal values of 
the computed solution (e.g., see [31]). A hierarchical basis error estimator can also be used to 
approximate the Hessian of the exact solution. It is shown in [34] that the least-squares fitting and 
the hierarchical basis methods work comparably for all considered cases except for one where the 
diffusion coefficient is discontinuous and the interfaces are predefined in the mesh. In this case, the 
latter works better than the former since hierarchical basis estimation does not over-concentrate 
mesh elements near the interfaces. Since our main goal is to study DMP satisfaction instead of 
the discontinuity of the diffusion coefficient, we choose to use the least squares fitting method for 
Hessian recovery in our computation due to its simplicity and problem independent feature. 



It is interesting to note that the term in the bracket in (60) can be viewed as a Riemann sum 
of an integral, i.e., 



E 

KeT h 



\K\B 



d+2 

K 



det (D)' 



i 

' d+2 



0-1 



d 

I d+2 



\H(u) 



2d 

d +2dx. 



Thus, the interpolation error has an asymptotic bound as 



\u - LJ,,«| H i ( n) < CN-~ d \ \K\B 



~ CN^s [ I det (D)~dT2 llp™ 1 !!^ 




\H(u)\\\*+idx 



d+2 
2d 



(61) 



We emphasize that both the satisfaction of DMP and mesh adaptation (through minimization 
of an error bound) are taken into account in the definition of metric tensor ( |55[ ). An interesting 
question is what the interpolation error bound looks like if mesh adaptation is not taken into 
consideration. For example, we consider a case 9k = 1 in (39). This gives the metric tensor 



M K 



(62) 



Recall that the interpolation error is bounded in (48), i.e. 



<C [ Y, \ k \~ b k 

,K&T h 



(63) 



where Bk is defined in (49). Moreover, for an M-uniform mesh corresponding to this metric tensor 

(64) 



the equidistribution condition ( 35 ) reduces to 



det(Bjf)-a|^T| 



N ' 
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where Oh= det(D^) 2 \K\. Inserting (64) into (63), we have 
KeT h 



V- h u\ mm < C \J2 \ K \ (det(B K )^y B K 
\KeT h 

= CN--da\ I \K\det(B K )^B K 



CN- L d det(B K )-^\K\ ^ \K\det(B K )*B K 



J<eT h 



.KeT h 



Thus, 



n h u\ H i {n) < CN- 1 * Yl det(® K )-*\K 




^ \K\det(B K )^B K 



'H ■ \\B\H(u)\\\ 2 dx 



(65) 
(66) 



This is the interpolation error bound for an M -uniform mesh corresponding to metric tensor ( 62 ) 
Prom Holder's inequality, it follows that 




< [ det(B K )~^\K\ ^ \K\det(B K )^B K 



Thus, the solution-dependent factor of bound (60) is small than or equal to that of bound (65). In 



this sense, MpMP+adap defined in (55) leads to a more accurate interpolant than Mdmp defined 
in d62]) (or ((39} with 9 K = 1). 



Moreover, from the standard interpolation theory we recall that the interpolation error for a 
uniform mesh is bounded by 



u ■ 



< CN~ 



\V 2 u\\ 2 dx 



(67) 



It is easy to see that the solution dependent factor in error bound (61) for MpMP+adap is in 
the order of |V 2 u 



2d and those in (66) for Mdmp and (67) for a uniform mesh are in the 

LJ+2 (fi) 



order of |V u|^2(q). Thus, (61) has the smallest solution dependent factor, an indication of the 



advantage of using adaptive meshes. On the other hand, the error bounds (61) and (66) depend 
on the determinant and norm of the diffusion matrix D and its inverse. This indicates that DMP 
satisfaction may sacrifice accuracy. Indeed, as we shall see in the next section, the solution error 
for DMP-bound meshes can sometimes be larger than that for a uniform mesh. 
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Given a mesh 




Solve PDE 




Recover solution 
derivatives and 
compute M 




Generate 
new mesh 
according to M 


i 


> 


> 








iteration 









Figure 2: An iterative procedure for adaptive mesh solution of PDEs 



5 Numerical results 



In this section we present three two-dimensional examples to demonstrate the performance of metric 



tensors Mdmp i n (39) with 6k = 1 based on DMP satisfaction and MpMP+adap m (55) combining 



DMP satisfaction and mesh adaptivity. For comparison purpose, we also include numerical results 
obtained with almost uniform meshes (labelled with M un if) and with a metric tensor M a d ap based 
on minimization of a bound on the H 1 semi-norm of linear interpolation error 



M, 



where 



adap,K 



PK 



p d K det 



I+—\H K (u) 



I+—\H K (u) 
ah 



(68) 



I+—\H K (u) 
ah 



d 
d + 2 



dct 



I+—\H K {u) 
ah 



i 

d+2 



and Qh is defined implicitly through 



\K\p K = 2|n| 

KeT h 



Once again, the second derivatives of the exact solution are replaced in actual computation with 
approximations obtained with a Hessian recovery technique (using piecewise quadratic polynomials 
fitting in least-squares sense to nodal values of the computed solution |31j). 

An iterative procedure for solving PDEs is shown in Fig. [2] In the current computation, each 
run is stopped after ten iterations. We have found that there is very little improvement in the 
computed solution after ten iterations for all the examples considered. A new mesh is generated 
using the computer code BAMG (bidimensional anisotropic mesh generator) developed by Hecht 
[29] based on a Delaunay-type triangulation method [TT]. The code allows the user to supply 
his/her own metric tensor defined on a background mesh. In our computation, the background 
mesh has been taken as the most recent mesh available. 



Example 5.1 The first example is to consider BVP ([I]) and ^ with 



/ = 0, ft=[0,l] 2 \ 



4 5 
9' 9 



9 = 0on T ou t, 9 = 2 on Tj 



where T out and Ti n are the outer and inner boundaries of Q, respectively; see Fig. [3| The diffusion 



matrix is given by (31) with k\ = 1000, &2 = 1, and 9 being the angle of the primary diffusion 



direction (parallel to the first eigenvector of . 
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This example satisfies the maximum principle and the solution (whose analytical expression is 
unavailable) stays between and 2. Our goal is to produce a numerical solution which also satisfies 
DMP and stays between and 2. Moreover, for both cases with a constant and a variable 9 we 
consider, the exact solution has sharp jumps near the inner boundary (cf. Figs. [4] and [8]) so mesh 
adaptation is needed for a proper resolution of them. This example has been studied in [33l [50] . 

We first consider the case of constant D with 9 = 7r/4. Fig. [4] shows finite element solutions 
obtained with M uni f and MpMP+adap- Meshes and solution contours obtained with various metric 
tensors are shown in Figs. [5]and[6j respectively. No overshoots in the finite element solutions are 
observed for all cases. However, undershoots and unphysical minima occur in the solutions obtained 
with M unif (u min = -0.0602) and M adap (u min = -0.0039) (cf. Fig. g(a) and (b)). Fig. [7] shows 
the decrease of —u m i n as the mesh is refined. For the range of the number of mesh elements 
considered, the undershooting improves at a rate of —u m i n = O(N~ ' 5 ) for both M un if and M a d ap . 
On the other hand, the results confirm the theoretical prediction that the solutions obtained with 
Mdmp and MpMP+adap satisfy DMP and no overshoot /undershoot and no unphysical extremum 
occur. It should be pointed out that the solution contour obtained with an almost uniform mesh is 
very smooth but the sharp jumps of the solution are smeared; see Figs. |4^a) and[6](a). The solution 
contours obtained with Mdmp and MpMP+adap are comparable to the one obtained with M adav . 

Next we consider a case of variable D with 9 = tt sin(x) cos(y). The finite element solutions, 
meshes, and solution contours are shown in Figs. [8j [9| and [TUJ respectively. Similar observations 
as for the constant D case can be made. Especially, undershoots and unphysical extrema occur in 
the solutions obtained with M un if and M a d ap but not with Mdmp and MpMP+adap- Once again, 
the results confirm our theoretical predictions in the previous sections. 



u 



= 2 



u = 



Figure 3: The physical domain and boundary conditions for Example 5.1 



Example 5.2 In this exampL 



.e, we consider BVP ^ and ^ with 



/ = 0, g(x,0)=g(16,y) 
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Figure 4: Example 

MoMP+adap- 



(a): M uni f (b): M DMP+adap 

5.1 with constant D. Finite element solutions obtained with (a) M un if and (b) 



The diffusion matrix is defined as 



D(a;,y) 



500.5 499.5 
499.5 500.5 



This is a simple example with a constant but anisotropic D and with a continuous boundary 
condition. It satisfies the maximum principle and its solution stays between and 1. 



Numerical solutions, meshes, and solution contours are shown in Figs. 11, 12, and 13, respec- 
tively. For this example, both undershoots and overshoots are observed in the computed solutions 
with M uni f and M ac [ ap but not with with Mdmp and MpMP+adap- This example demonstrates that 
a scheme violating DMP can produce unphysical extrema even for a simple problem with constant 
diffusion, continuous boundary conditions, and a convex domain. 



Example 5.3 This example is given by @ and @ with 

fi = (0,1) x (0,1), f(x,y) = | -f ^ ^ q 5 ' u = u exact on dtt, 

B( X ,y) = l Du ifx< °- 5 ' Dl =( 1 °V ^=f 10 3 V 
v ,y) \ d 2 , if x > o.5, yo 1 J' v 3 1 y 

The problem has the exact solution 

. / 1 - 2y 2 + Axy + 2y + 6x, if x < 0.5 

u(x,y) = < „ (d9 
v y | -2y 2 + 1.6a;y-0.6x + 3.2y + 4.3, if x > 0.5. v ; 

Note that the value and primary diffusion direction of the diffusion matrix change across the line 
x = 0.5. This example has been studied in [33]. 
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M, 



unif j 



N = 2460 



(b): M adap , N = 2583 





(c): M DM p, N = 2530 
Figure 5: Example |5 . 1 1 with constant 



(d): M DMP+adap , N = 2474 
Meshes obtained with different metric tensors. 



Solutions and meshes obtained with various metric tensors are shown in Fig. 14 For this 
example, no overshoots and undershoots are observed for all numerical solutions. The meshes 
obtained with Mdmp and MpMP+adap show a better alignment with the primary diffusion direction 
than that obtained with M ada p- Moreover, elements are concentrated along the line x = 0.5 for the 
meshes obtained with M ada p and M^MP+adap whereas there is no concentration in the mesh shown 



in Fig. 14 d) for Mdmp- The results are consistent with what is expected from the construction 



of the metric tensors. 

The exact solution is available for this example. The H l semi-norm and L 2 norm of the error 



are shown in Fig. 15 as functions of the number of mesh elements. Metric tensor M ac i a p leads to 
far more accurate results than the other three metric tensors, which produce comparable results 
for the considered range of N. Moreover, M adav and MpMP+adap gi ye the same convergence 
rate, i.e., \e h \ m{n) = 0(N-°- 5 ) and \\e h \\ L 2 (n: 
slower convergence rate, |e"|#i(Q) = 0(iV -a25 ) and 
advantage of using adaptive meshes. Interestingly, the results in [13] (Table 4) obtained with a 



0(N ), while M uni f and Mdmp result in a 
\l 2 (q) = 0(-/V~°- 5 ). This demonstrates the 
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(a): M unif , Umin = -0.0602 (b): M adap , u min = -0.0039 




(c): Mdmp, Umin = (d): M DMP+adap , u min = 

Figure 6: Example |5.1| with constant D. Contours of the finite element solutions obtained with 
different metric tensors. 



slope-limited scheme for triangular meshes also show a similar slow convergence. 

It should be pointed out that the above results have been obtained when the interface (x = 0.5) 
is not predefined in the mesh. If the interface is predefined in the mesh, then the solution ( |69| ) 
can be approximated accurately in the linear finite element space. As shown in Fig. \TE[ all metric 
tensors produce comparable solutions and the same convergence rate |e fe |^i(Q) = O(N~ ' 5 ) and 
\\e h \\Li(n) = 0(N~ 1 ). 



6 Conclusions and comments 

In the previous sections we have developed a mesh condition ( p4| ) under which the linear finite ele- 
ment approximation of anisotropic diffusion problem ([!]) and ^ validates the discrete counterpart 
of the maximum principle satisfied by the continuous problem. The condition is a generalization of 
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0.1 



000 10000 100000 

N: number o( elements 



Figure 7: Example 5.1 with constant D. The undershoot, —u m i n , is shown as functions of the 
number of elements. 





fa): M 



unif 



(b): M^MP+adap 



Figure 8: Example 5.1 with variable D. Finite element solutions obtained with (a) M un if and (b) 

MoMP+adap- 



the well known non-obtuse angle condition developed for isotropic diffusion problems and requires 
that the dihedral angles of mesh elements measured in a metric depending only on the diffusion 
matrix be non-obtuse. 



We have also developed two variants of the anisotropic non-obtuse angle condition, (33) and 



(34), which can be more convenient to use in actual mesh generation. Indeed, metric tensor (39) 



for use in anisotropic mesh generation is derived based on ( 34 ) for accounting for DMP satisfaction 



Moreover, an optimal metric tensor ( 55 ) accounting for both DMP satisfaction and mesh adaptation 



is obtained from (34) by minimizing an interpolation error bound. Features of these metric tensors 



are illustrated in numerical examples. 



It is worth pointing out that condition (24) has been derived based on the local stiffness matrix 



on a mesh element. Like the non-obtuse angle condition for isotropic diffusion problems, (24) 



may be relaxed by considering the global stiffness matrix as a whole [19]. Moreover, we have 
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(a): M umf , iV = 2460 



(b): M adap , N = 2568 




(c): M DM p, N = 2510 



Figure 9: Example 5.1 with variable 



(d): M DMP+adap , N = 2463 
Meshes obtained from different metric tensors. 



restricted our attention to linear PDE ([I]) and Dirichlet boundary condition ([2]). But the procedure 
developed in this work can be extended to problems with nonlinear diffusion D = D(x,u, Vu) and 
mixed boundary conditions (e.g., see [3^1 E3 EHl IB] ) without major modification. 

Although the numerical examples have been presented in 2D, the anisotropic non-obtuse angle 



condition ([24]) and the corresponding metric tensor formulas (39), (55), and (62) are d-dimensional 
(d = 1,2,3). In 3D, a Delaunay triangulation may not guarantee the satisfaction of DMP [49|. 
Nevertheless, some polyhedrons can be decomposed into tetrahedra satisfying the non-obtuse angle 
condition (30) and therefore the numerical solution satisfies DMP; e.g., see [33]. It is expected that 



this will also work for the anisotropic non-obtuse angle condition (24) for a given metric tensor M. 



On the other hand, the existence of the decomposition of an arbitrary polyhedron into non-obtuse 
tetrahedra is an open problem [33] . It is also unclear if a 3D triangulation can be generated to 



(approximately) satisfy the M-uniform mesh conditions ( 35 and 36 ) . Those are interesting topics 
to investigate in the future. 
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(c): Mdmp, Umin = (d): M DMP+adap , u m i n = 

Figure 10: Example |5.1| with variable D. Contours of the finite element solutions obtained with 
different metric tensors. 



Acknowledgment. The work was supported in part by the National Science Foundation 
(USA) under grant DMS-0712935. 

References 

[1] I. Aavatsmark, T. Barkve, 0. B0e, and T. Mannseth. Discretization on unstructured grids 
for inhomogeneous, anisotropic media. I. Derivation of the methods. SI AM J. Sci. Comput., 
19:1700-1716 (electronic), 1998. 

[2] I. Aavatsmark, T. Barkve, 0. B0e, and T. Mannseth. Discretization on unstructured grids 
for inhomogeneous, anisotropic media. II. Discussion and numerical results. SIAM J. Sci. 
Comput, 19:1717-1736 (electronic), 1998. 



25 



(a): M adap , u min = -0.0195 (b): M DMP+adap , u min = 

Fi gure 11; Example |5.2| Finite element solutions obtained with (a) M ada p and (b) M-DMP+adap' 



[3] D. Ait-Ali-Yahia, G. Baruzzi, W. G. Habashi, M. Fortin, J. Dompierre, and M.-G. Val- 
let. Anisotropic mesh adaptation: towards user- independent, mesh-independent and solver- 
independent CFD. Part II: Structured grids. Int. J. Numer. Meth. Fluids, 39:657-673, 2002. 

[4] H. Borouchaki, P. L. George, P. Hecht, P. Laug, and E. Saletl. Delaunay mesh generation 
governed by metric specification: Part I. Algorithms. Fin. Elem. Anal. Des., 25:61-83, 1997. 

[5] H. Borouchaki, P. L. George, and B. Mohammadi. Delaunay mesh generation governed by 
metric specification: Part II. Applications. Fin. Elem. Anal. Des., 25:85-109, 1997. 

[6] F. J. Bossen and P. S. Heckbert. A pliant method for anisotropic mesh generation. In Pro- 
ceedings, 5th International Meshing Roundtable, pages 63-74, Sandia National Laboratories, 
Albuquerque, NM, 1996. Sandia Report 96-2301. 

[7] J. U. Brackbill and J. S. Saltzman. Adaptive zoning for singular problems in two dimensions. 
J. Comput. Phys., 46:342-368, 1982. 

[8] J. Brandts, S. Korotov, and M. Kffzek. Dissection of the path-simplex in W 1 into n path- 
subsimplices. Lin. Alg. AppL, 421:382-393, 2007. 

[9] J. Brandts, S. Korotov, and M. Kffzek. The discrete maximum principle for linear simplicial 
finite element approximations of a reaction-diffusion problem. Lin. Alg. AppL, 429:2344-2357, 
2008. 

[10] E. Burman and A. Ern. Discrete maximum principle for Galerkin approximations of the 
Laplace operator on arbitrary meshes. C. R. Acad. Sci. Paris, Ser.I 338:641-646, 2004. 

[11] M. J. Castro-Diaz, F. Hecht, B. Mohammadi, and O. Pironneau. Anisotropic unstructured 
mesh adaption for flow simulations. Int. J. Numer. Meth. Fluids, 25:475-491, 1997. 



26 



(a): M unif , iV = 2 480 



(b): M adap , N = 2480 




2 4 6 8 10 12 14 16 




2 4 6 8 10 12 14 16 



(c): MflMP, N = 2754 (d): M DMP+adap , N = 2529 

Figure 12: Example |5.2| The adaptive meshes obtained with various metric tensors. 

[12] T. F. Chan and J. Shen. Non-texture inpainting by curvature driven diffusions (CDD). J. Vis. 
Commun. Image Rep, 12:436-449, 2000. 

[13] T. F. Chan, J. Shen, and L. Vese. Variational PDE models in image processing. Not. AMS 
J., 50:14-26, 2003. 

[14] P. G. Ciarlet. Discrete maximum principle for finite difference operators. Aequationes Math., 
4:338-352, 1970. 

[15] P. G. Ciarlet and P.-A. Raviart. Maximum principle and uniform convergence for the finite 
element method. Comput. Meth. Appl. Mech. Eng., 2:17-31, 1973. 

[16] P. I. Crumpton, G. J. Shaw, and A. F. Ware. Discretisation and multigrid solution of elliptic 
equations with mixed derivative terms and strongly discontinuous coefficients. J. Comput. 
Phys., 116:343-358, 1995. 



27 



(a): M uni f, u min = -0.0280, u max = 1.0209 (b): M adap , u min = -0.0195, u max = 1.0160 




(c): MdMP, Umin = 0, U m ax = (d): M DMP+adap , U m in = 0, U m ax = 

Figure 13: Example |5.2| Contours of the finite element solutions obtained with different metric 
tensors. 

[17] J. Dompierre, M.-G. Vallet, Y. Bourgault, M. Fortin, and W. G. Habashi. Anisotropic mesh 
adaptation: towards user-independent, mesh- independent and solver-independent CFD. Part 
III: Unstructured meshes. Int. J. Numer. Meth. Fluids, 39:675-702, 2002. 

[18] A. Draganescu, T. F. Dupont, and L. R. Scott. Failure of the discrete maximum principle for 
an elliptic finite element problem. Math. Comp., 74:1-23, 2004. 

[19] A. S. Dvinsky. Adaptive grid generation from harmonic maps on riemannian manifolds. J. 
Comput. Phys., 95:450-476, 1991. 

[20] G. T. Eigestad, I. Aavatsmark, and M. Espedal. Symmetry and M-matrix issues for the 
O-method on an unstructured grid. Comput. Geosci., 6:381-404, 2002. 

[21] A. Ern and J. L. Guermond. Theory and Practice of Finite Elements. Sprigner-Verlag, New 
York, 2004. 



28 



[22] T. Ertekin, J. H. Abou-Kassem, and G. R. King. Basic Applied Reservoir Simulation. SPE 
textbook series, Vol. 7, Richardson, Texas, 2001. 

[23] P. Frey and P. L. George. Mesh Generation: Application to Finite Elements. Hermes Science, 
Oxford and Paris, 2000. 

[24] R. V. Garimella and M. S. Shephard. Boundary layer meshing for viscous flows in complex 
domain. In Proceedings, 7th International Meshing Roundtable, pages 107-118, Sandia National 
Laboratories, Albuquerque, NM, 1998. 

[25] S. Gunter and K. Lackner. A mixed implicit-explicit finite difference scheme for heat transport 
in magnetised plasmas. J. Comput. Phys., 228:282-293, 2009. 

[26] S. Gunter, K. Lackner, and C. Tichmann. Finite element and higher order difference formula- 
tions for modelling heat transport in magnetised plasmas. J. Comput. Phys., 226:2306-2316, 
2007. 

[27] S. Gunter, Q. Yu, J. Kruger, and K. Lackner. Modelling of heat transport in magnetised 
plasmas using non-aligned coordinates. J. Comput. Phys., 209:354-370, 2005. 

[28] W. G. Habashi, J. Dompierre, Y. Bourgault, D. Ait-Ali-Yahia, M. Fortin, and M.-G. Val- 
let. Anisotropic mesh adaptation: towards user-independent, mesh-independent and solver- 
independent CFD. Part I: General principles. Int. J. Numer. Meth. Fluids, 32:725-744, 2000. 

[29] F. Hecht. Bidimensional anisotropic mesh generator. Technical report, INRIA, Rocquencourt, 
1997. Source code: http://www.ann.jussieu.fr/ hecht /ftp/bamg/. 

[30] W. Huang. Variational mesh adaptation: isotropy and equidistribution. J. Comput. Phys., 
174:903-924, 2001. 

[31] W. Huang. Metric tensors for anisotropic mesh generation. J. Comput. Phys., 204:633-665, 
2005. 

[32] W. Huang. Mathematical principles of anisotropic mesh adaptation. Comm. Comput. Phys., 
1:276-310, 2006. 

[33] W. Huang and X. P. Li. An anisotropic mesh adaptation method for the finite element solution 
of variational problems. Fin. Elem. Anal. Des., 46:61-73, 2010. 

[34] W. Huang and L. Kamenski and J. Lang. A new anisotropic mesh adaptation method based 
upon hierarchical a posteriori error estimates. J. Comput. Phys., 229:2179-2198, 2010. 

[35] Olivier-P. Jacquotte. A mechanical model for a new grid generation method in computational 
fluid dynamics. Comput. Meth. Appl. Mech. Engrg., 66:323-338, 1988. 

[36] J. Karatson and S. Korotov. Discrete maximum principles for finite element solutions of 
nonlinear elliptic problems with mixed boundary conditions. Numer. Math., 99:669-698, 2005. 

[37] J. Karatson and S. Korotov. Discrete maximum principles for finite element solutions of some 
mixed nonlinear elliptic problems using quadratures. J. Comput. Appl. Math., 192:75-88, 2006. 



29 



[38] J. Karatson, S. Korotov, and M. Kffzek. On discrete maximum principles for nonlinear elliptic 
problems. Math. Comput. Sim., 76:99-108, 2007. 

[39] D. A. Karras and G. B. Mertzios. New PDE-based methods for image enhancement using SOM 
and Bayesian inference in various discretization schemes. Meas. Sci. TechnoL, 20:104012, 2009. 

[40] P. Knupp, L. Margolin, and M. Shashkov. Reference jacobian optimization-based rezone strate- 
gies for arbitrary lagrangian eulerian methods. J. Comput. Phys., 176:93-128, 2002. 

[41] P. M. Knupp. Jacobian-weighted elliptic grid generation. SI AM J. Sci. Comput., 17:1475-1490, 
1996. 

[42] R. Kornhuber and R. Roitzsch. On adaptive grid refinement in the presence of internal or 
boundary layers. IMPACT Comput. Sci. Engrg., 2:40-72, 1990. 

[43] D. Kuzmin, M. J. Shashkov, and D. Svyatskiy. A constrained finite element method satisfying 
the discrete maximum principle for anisotropic diffusion problems. J. Comput. Phys., 228:3448- 
3463, 2009. 

[44] M. Kffzek and Q. Lin. On diagonal dominance of stiffness matrices in 3D. East-West J. 
Numer. Math., 3:59-69, 1995. 

[45] J. Lang. An adaptive finite element method for convection-diffusion problems by interpolation 
techniques. Technical Report TR 91-4, Konrad-Zuse-Zentrum Berlin, 1991. 

[46] C. Le Potier. Schema volumes finis monotone pour des operateurs de diffusion fortement 
anisotropes sur des maillages de triangles non structures. C. R. Math. Acad. Sci. Paris, 
341:787-792, 2005. 

[47] C. Le Potier. A nonlinear finite volume scheme satisfying maximum and minimum principles 
for diffusion operators. Int. J. Finite Vol., 6:20, 2009. 

[48] C. Le Potier. Un schema lineaire verifiant le principe du maximum pour des operateurs de 
diffusion tres anisotropes sur des maillages deformes. C. R. Math. Acad. Sci. Paris, 347:105- 
110, 2009. 

[49] F. W. Letniowski. Three-dimensional delaunay triangulations for finite element approximations 
to a second-order diffusion operator. SIAM J. Sci. Stat. Comput., 13:765-770, 1992. 

[50] X. P. Li, D. Svyatskiy, and M. Shashkov. Mesh adaptation and discrete maximum principle 
for 2D anisotropic diffusion problems. Technical report, LANL, 2007. Final Report of the 
Summer Student Program. 

[51] K. Lipnikov, M. Shashkov, D. Svyatskiy, and Yu. Vassilevski. Monotone finite volume schemes 
for diffusion equations on unstructured triangular and shape-regular polygonal meshes. J. 
Comput. Phys., 227:492-512, 2007. 

[52] R. Liska and M. Shashkov. Enforcing the discrete maximum principle for linear finite element 
solutions of second-order elliptic problems. Comm. Comput. Phys., 3:852-877, 2008. 



30 



[53] M. J. Mlacnik and L. J. Durlofsky. Unstructured grid optimization for improved monotonicity 
of discrete solutions of elliptic equations with highly anisotropic coefficients. J. Comput. Phys., 
216:337-361, 2006. 

[54] D. Mumford and J. Shah. Optimal approximations by piecewise smooth functions and associ- 
ated variational problems. Commun. Pure Appl. Math, 42:577-685, 1989. 

[55] K. Nishikawa and M. Wakatani. Plasma Physics. Springer- Verlag Berlin Heidelberg, New 
York, 2000. 

[56] J. Peraire, M. Vahdati, K. Morgan, and O. C. Zienkiewicz. Adaptive remeshing for compressible 
flow computations. J. Comput. Phys., 72:449-466, 1987. 

[57] P. Perona and J. Malik. Scale-space and edge detection using anisotropic diffusion. IEEE 
Trans. Pattern Anal. Mach. Intel., 12:629-639, 1990. 

[58] W. Rachowicz. An anisotropic /i-type mesh refinement strategy. Comput. Meth. Appl. Mech. 
Engrg., 109:169-181, 1993. 

[59] W. Rachowicz. An anisotropic /i-adaptive finite element method for compressible Navier-Stokes 
equations. Comput. Meth. Appl. Mech. Engrg., 146:231-252, 1997. 

[60] J. Remacle, X. Li, M. S. Shephard, and J. E. Flaherty. Anisotropic adaptive simulation of 
transient flows using discontinuous Galerkin methods. Int. J. Numer. Meth. Engr., 62:899-923, 
2003. 

[61] P. Sharma and G.W. Hammett. Preserving monotonicity in anisotropic diffusion. J. Comput. 
Phys., 227:123-142, 2007. 

[62] D. M. Y. Sommerville. An Introduction to the Geometry of n Dimensions. Methuen & Co. 
LTD., London, 1929. 

[63] T.H. Stix. Waves in Plasmas. Amer. Inst. Phys., New York, 1992. 

[64] G. Stoyan. On a maximum principle for matrices, and on conservation of monotonicity. With 
applications to discretization methods. Z. Angew. Math. Mech., 62:375-381, 1982. 

[65] G. Stoyan. On maximum principles for monotone matrices. Lin. Alg. Appl., 78:147-161, 1986. 

[66] G. Strang and G. J. Fix. An Analysis of the Finite Element Method. Prentice Hall, Englewood 
Cliffs, NJ, 1973. 

[67] R. S. Varga. Matrix Iterative Analysis. Prentice-Hall, New Jersey, 1962. 

[68] R. S. Varga. On a discrete maximum principle. SIAM J. Numer. Anal, 3:355-359, 1966. 

[69] J. Weickert. Anisotropic Diffusion in Image Processing. Teubner- Verlag, Stuttgart, Germany, 
1998. 

[70] A. M. Winslow. Adaptive mesh zoning by the equipotential method. Technical Report UCID- 
19062, Lawrence Livemore Laboratory, 1981. (unpublished). 



31 



[71] J. Xu and L. Zikatanov. A monotone finite element scheme for convection-diffusion equations. 
Math. Comput., 69:1429-1446, 1999. 

[72] S. Yamakawa and K. Shimada. High quality anisotropic tetrahedral mesh generation via 
ellipsoidal bubble packing. In Proceedings, 9th International Meshing Roundtable, pages 263- 
273, Sandia National Laboratories, Albuquerque, NM, 2000. Sandia Report 2000-2207. 



32 





(a): M a d ap , numerical solution, u m i n = 



(b): M adap , mesh, iV = 2362 



Jit* 
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(c): MdmPi numerical solution, u m i n = 



(d): M DM p, mesh, N = 2415 





(e): M DM p +adapt numerical solution, u m i n — (f): M DM p+adap, mesh, N = 2490 

Figure 14: Example 5.3 Numerical solutions and meshes obtained with three metric tensors. 
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Figure 15: Example 5.3. The H 1 semi-norm and L 2 norm of solution error are shown as functions 
of the number of elements for metric tensors M un if, M a d a p, MdatPj and MjyMP+adap- 
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Figure 16: Example 5.3 The H l semi-norm and L 2 norm of solution error are shown as functions 
of the number of elements for metric tensors M un if, M a d api Mdmp, and MoMP+adap- The interface 
(x = 0.5) is predefined in the mesh. 
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